# Load libraries, install if needed
pkgs <- c("tidyverse", "tidylog", "RCurl", "devtools", "modelr", "ggtext",
"viridis", "RColorBrewer", "here", "sdmTMBextra")
if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){
install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
}
invisible(lapply(pkgs, library, character.only = T))
# Packages not on CRAN or dev version
# remotes::install_github("pbs-assess/sdmTMB", dependencies = TRUE)
library(sdmTMB)
# Source code for map plots
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/R/functions/map-plot.R")
options(ggplot2.continuous.colour = "viridis")
# Set path
home <- here::here()Fit tweedie/delta models to biomass density
Intro
Fit climate-agnostic sdms for calculating biomass trends and velocities
Load packages & source functions
Read data
# Read data
d <- #readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/catch_clean.csv") %>%
readr::read_csv(paste0(home, "/data/clean/catch_clean.csv")) %>%
rename(X = x, Y = y) %>%
pivot_longer(c(cod_adult, cod_juvenile, dab_adult, dab_juvenile, flounder_adult,
flounder_juvenile, plaice_adult, plaice_juvenile),
names_to = "group", values_to = "density") %>%
separate(group, into = c("species", "life_stage"), remove = FALSE) %>%
drop_na(depth, temp, oxy, sal, density)Rows: 12385 Columns: 25
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (4): country, haul_id, ices_rect, month_year
dbl (21): year, lat, lon, quarter, month, sub_div, oxy, temp, sal, depth, x,...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed 2 variables (X, Y)
pivot_longer: reorganized (cod_adult, cod_juvenile, dab_adult, dab_juvenile, plaice_adult, …) into (group, density) [was 12385x25, now 99080x19]
drop_na: removed 9,570 rows (10%), 89,510 rows remainingd <- d %>% filter(!species == "dab")filter: removed 23,124 rows (26%), 66,386 rows remainingScale variables
d <- d %>%
mutate(temp_sc = as.vector(scale(temp)),
temp_sq = temp_sc^2,
oxy_sc = as.vector(scale(oxy)),
oxy_sq = oxy_sc^2,
sal_sc = as.vector(scale(sal)),
depth_sc = as.vector(scale(depth)),
depth_sq = depth_sc*depth_sc,
quarter_f = as.factor(quarter),
year_f = as.factor(year),
year_ct = year - mean(year),
year_sc = (year - mean(year)) / sd(year)) %>%
filter(quarter == 4)mutate: new variable 'temp_sc' (double) with 12,194 unique values and 0% NA new variable 'temp_sq' (double) with 12,194 unique values and 0% NA new variable 'oxy_sc' (double) with 12,192 unique values and 0% NA new variable 'oxy_sq' (double) with 12,192 unique values and 0% NA new variable 'sal_sc' (double) with 12,195 unique values and 0% NA new variable 'depth_sc' (double) with 4,519 unique values and 0% NA new variable 'depth_sq' (double) with 4,519 unique values and 0% NA new variable 'quarter_f' (factor) with 2 unique values and 0% NA new variable 'year_f' (factor) with 29 unique values and 0% NA new variable 'year_ct' (double) with 29 unique values and 0% NA new variable 'year_sc' (double) with 29 unique values and 0% NAfilter: removed 40,246 rows (61%), 26,140 rows remainingRead and scale the prediction grid
pred_grid <- bind_rows(readr::read_csv(paste0(home, "/data/clean/pred_grid_(1_2).csv")),
readr::read_csv(paste0(home, "/data/clean/pred_grid_(2_2).csv")))Rows: 457436 Columns: 14
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): month_year, ices_rect
dbl (12): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_div
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 490110 Columns: 14
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): month_year, ices_rect
dbl (12): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_div
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# Scale variables with respect to data, not the prediction grid!
pred_grid <- pred_grid %>%
drop_na(depth, temp, oxy, sal) %>%
mutate(temp_sc = (temp - mean(d$temp)) / sd(d$temp),
temp_sq = temp_sc^2,
oxy_sc = (oxy - mean(d$oxy)) / sd(d$oxy),
oxy_sq = oxy_sc^2,
sal_sc = (sal - mean(d$sal)) / sd(d$sal),
depth_sc = (depth - mean(d$depth)) / sd(d$depth),
depth_sq = depth_sc*depth_sc,
quarter_f = as.factor(quarter),
year_ct = year - mean(d$year),
year_sc = (year - mean(d$year)) / sd(d$year)) %>%
filter(quarter == 4)drop_na: removed 4,118 rows (<1%), 943,428 rows remaining
mutate: new variable 'temp_sc' (double) with 917,493 unique values and 0% NA
new variable 'temp_sq' (double) with 917,493 unique values and 0% NA
new variable 'oxy_sc' (double) with 914,731 unique values and 0% NA
new variable 'oxy_sq' (double) with 914,731 unique values and 0% NA
new variable 'sal_sc' (double) with 897,794 unique values and 0% NA
new variable 'depth_sc' (double) with 7,572 unique values and 0% NA
new variable 'depth_sq' (double) with 7,572 unique values and 0% NA
new variable 'quarter_f' (factor) with 2 unique values and 0% NA
new variable 'year_ct' (double) with 29 unique values and 0% NA
new variable 'year_sc' (double) with 29 unique values and 0% NA
filter: removed 471,714 rows (50%), 471,714 rows remainingFit the model that is preferred by most (m3, i.e., temp + temp_sq + bp_oxyg)
Use only quarter 4 since that’s what we are looking for (warmer and less oxygen than q1), and not also we use a spatial trend model Plot residuals, save model object, plot conditional effects!
# https://pbs-assess.github.io/sdmTMB/articles/spatial-trend-models.html
pred_grid_list <- list()
res_list <- list()
for(i in unique(d$group)) {
dd <- d %>%
filter(group == i & quarter == 4)
mesh <- make_mesh(dd, xy_cols = c("X", "Y"), cutoff = 15)
plot(mesh)
# This is a variant of the main model without environmental fixed effects
# FIXME: finalize model structure poisson-link...
m <- sdmTMB(density ~ 0 + year_sc + depth_sc + depth_sq,
data = dd,
mesh = mesh,
family = delta_lognormal(type = "poisson-link"), # because we want to combine the spatial trends!
spatiotemporal = "IID",
spatial = "on",
spatial_varying = ~0 + year_sc,
time = "year")
# Check the model converges
sanity(m)
tidy(m, "ran_pars", conf.int = TRUE, model = 1)
#tidy(m, "ran_pars", conf.int = TRUE, model = 2)
# Plot residuals
# dd$res1 <- residuals(m_test3, model = 1)
# dd$res2 <- residuals(m_test2, model = 2)
#
# qqnorm(dd$res1); qqline(dd$res1)
# qqnorm(dd$res2[is.finite(dd$res2)]); qqline(dd$res2[is.finite(dd$res2)])
# MCMC
samps1 <- sdmTMBextra::predict_mle_mcmc(m, mcmc_iter = 201, mcmc_warmup = 200, model = 1)
samps2 <- sdmTMBextra::predict_mle_mcmc(m, mcmc_iter = 201, mcmc_warmup = 200, model = 2)
dd$mcmc_res1 <- as.vector(residuals(m, type = "mle-mcmc", mcmc_samples = samps1, model = 1))
dd$mcmc_res2 <- as.vector(residuals(m, type = "mle-mcmc", mcmc_samples = samps2, model = 2))
qqnorm(dd$mcmc_res1);qqline(dd$mcmc_res1)
qqnorm(dd$mcmc_res2[is.finite(dd$mcmc_res2)]); qqline(dd$mcmc_res2[is.finite(dd$mcmc_res2)])
# Store residuals
res_list[[i]] <- dd
# Predict on grid
p <- predict(m, newdata = pred_grid) %>%
mutate(group = i)
names(p)
plot_map_fc +
geom_raster(data = p, aes(X*1000, Y*1000, fill = zeta_s_year_sc1)) +
scale_fill_gradient2()
plot_map_fc +
geom_raster(data = p, aes(X*1000, Y*1000, fill = zeta_s_year_sc2)) +
scale_fill_gradient2()
# Save
pred_grid_list[[i]] <- p
}filter: removed 21,377 rows (82%), 4,763 rows remaining✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameters don't look unreasonably large
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.077508 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 775.08 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 7571.28 seconds (Warm-up)
Chain 1: 50.626 seconds (Sampling)
Chain 1: 7621.91 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.077073 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 770.73 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 7898.36 seconds (Warm-up)
Chain 1: 25.242 seconds (Sampling)
Chain 1: 7923.6 seconds (Total)
Chain 1: 
mutate: new variable 'group' (character) with one unique value and 0% NAfilter: removed 21,375 rows (82%), 4,765 rows remaining
✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameters don't look unreasonably large
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.076537 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 765.37 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 4969.02 seconds (Warm-up)
Chain 1: 12.507 seconds (Sampling)
Chain 1: 4981.52 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.078383 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 783.83 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 4781.52 seconds (Warm-up)
Chain 1: 38.625 seconds (Sampling)
Chain 1: 4820.15 seconds (Total)
Chain 1: 
mutate: new variable 'group' (character) with one unique value and 0% NAfilter: removed 22,374 rows (86%), 3,766 rows remaining
Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient
= tmb_obj$gr, : NA/NaN function evaluation✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameters don't look unreasonably large
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.063831 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 638.31 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 4166.92 seconds (Warm-up)
Chain 1: 19.534 seconds (Sampling)
Chain 1: 4186.45 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.061829 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 618.29 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 4168.19 seconds (Warm-up)
Chain 1: 19.385 seconds (Sampling)
Chain 1: 4187.58 seconds (Total)
Chain 1: 
mutate: new variable 'group' (character) with one unique value and 0% NAfilter: removed 22,373 rows (86%), 3,767 rows remaining
✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameters don't look unreasonably large
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.06067 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 606.7 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 5448.28 seconds (Warm-up)
Chain 1: 20.99 seconds (Sampling)
Chain 1: 5469.27 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.064953 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 649.53 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 5854.99 seconds (Warm-up)
Chain 1: 21.193 seconds (Sampling)
Chain 1: 5876.18 seconds (Total)
Chain 1: 
mutate: new variable 'group' (character) with one unique value and 0% NAfilter: removed 21,602 rows (83%), 4,538 rows remaining
Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient
= tmb_obj$gr, : NA/NaN function evaluation✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameters don't look unreasonably large
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.084936 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 849.36 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 7757.89 seconds (Warm-up)
Chain 1: 25.68 seconds (Sampling)
Chain 1: 7783.57 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.077033 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 770.33 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 8032.44 seconds (Warm-up)
Chain 1: 28.351 seconds (Sampling)
Chain 1: 8060.8 seconds (Total)
Chain 1: 
mutate: new variable 'group' (character) with one unique value and 0% NAfilter: removed 21,599 rows (83%), 4,541 rows remaining
✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameters don't look unreasonably large
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.083091 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 830.91 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 9653.7 seconds (Warm-up)
Chain 1: 51.438 seconds (Sampling)
Chain 1: 9705.14 seconds (Total)
Chain 1: Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceededWarning: Examine the pairs() plot to diagnose sampling problems
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.079155 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 791.55 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 9830.03 seconds (Warm-up)
Chain 1: 52.854 seconds (Sampling)
Chain 1: 9882.89 seconds (Total)
Chain 1: Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
Warning: Examine the pairs() plot to diagnose sampling problems
mutate: new variable 'group' (character) with one unique value and 0% NA
# Save predictions and sims as data frames
preds_grid <- dplyr::bind_rows(pred_grid_list) %>% as.data.frame()
res <- dplyr::bind_rows(res_list) %>% as.data.frame()Save
write_csv(res, paste0(home, "/output/residuals_04_sdms.csv"))
write_csv(preds_grid, paste0(home, "/output/data_pred_grids_04_random_sdms.csv"))Plot residuals
# Plot residuals
# FIXME: residuals look pretty crap - delta model instead?
res %>%
mutate(group2 = str_replace(group, "_", " "),
group2 = str_to_title(group2)) %>%
ggplot(aes(sample = mcmc_res1)) +
stat_qq(size = 0.75, shape = 21, fill = NA) +
facet_wrap(~group2) +
stat_qq_line() +
labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
theme(aspect.ratio = 1)mutate: new variable 'group2' (character) with 6 unique values and 0% NA
ggsave(paste0(home, "/figures/supp/qq_sdm_04_model1.pdf"), width = 17, height = 17, units = "cm")
res %>%
mutate(group2 = str_replace(group, "_", " "),
group2 = str_to_title(group2)) %>%
ggplot(aes(sample = mcmc_res2)) +
stat_qq(size = 0.75, shape = 21, fill = NA) +
facet_wrap(~group2) +
stat_qq_line() +
labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
theme(aspect.ratio = 1)mutate: new variable 'group2' (character) with 6 unique values and 0% NAWarning: Removed 12298 rows containing non-finite values (`stat_qq()`).Warning: Removed 12298 rows containing non-finite values (`stat_qq_line()`).
ggsave(paste0(home, "/figures/supp/qq_sdm_04_model2.pdf"), width = 17, height = 17, units = "cm")Warning: Removed 12298 rows containing non-finite values (`stat_qq()`).
Removed 12298 rows containing non-finite values (`stat_qq_line()`).